Designing 1D correlated-electron states by non-Euclidean topography of 2D monolayers

Two-dimensional (2D) bilayers, twisted to particular angles to display electronic flat bands, are being extensively explored for physics of strongly correlated 2D systems. However, the similar rich physics of one-dimensional (1D) strongly correlated systems remains elusive as it is largely inaccessible by twists. Here, a distinctive way to create 1D flat bands is proposed, by either stamping or growing a 2D monolayer on a non-Euclidean topography-patterned surface. Using boron nitride (hBN) as an example, our analysis employing elastic plate theory, density-functional and coarse-grained tight-binding method reveals that hBN’s bi-periodic sinusoidal deformation creates pseudo- electric and magnetic fields with unexpected spatial dependence. A combination of these fields leads to anisotropic confinement and 1D flat bands. Moreover, changing the periodic undulations can tune the bandwidth, to drive the system to different strongly correlated regimes such as density waves, Luttinger liquid, and Mott insulator. The 1D nature of these states differs from those obtained in twisted materials and can be exploited to study the exciting physics of 1D quantum systems.

R ecently, twisted bilayer graphene (TBG) at magic angles 1,2 and other van der Waals (vdW) heterostructures at small twist 3,4 have garnered great attention as material platforms for realizing 2D correlated physics with an unprecedented level of control. Several interesting electronic phases have been observed in these systems, such as correlated insulator 2,4 , superconductivity 1,5 , non-trivial electronic topology 6 , and magnetism 7,8 . Physically, these emergent phases can be attributed to the existence of the Bloch flat bands 9,10 , where the kinetic energy scale is quenched, and the role of electronic interactions is enhanced. The flat bands originate from the perturbation of the electronic structure by the longwavelength superlattice (moiré) period, which suppresses the group velocity in TBG at magic angles 11,12 and creates electronic confinement in other vdW heterostructures [13][14][15] . The moiré periods arise from either lattice mismatch or rotational misalignment between the layers with fine-tuning of the twist, posing challenges 16 in fabrication, variability between devices, and scalability. Moreover, examining the rich unexplored physics of 1D strongly correlated systems, likewise, is largely inaccessible by twists.
Below we describe creating flat bands through an alternative route, not requiring a twist angle. The strategy involves either growing 17 or stamping 18 a 2D material on a topographically patterned substrate with non-zero Gaussian curvature, that is non-Euclidean surface, Fig. 1a. To conform to such surfaces, a planar 2D crystal must deform, so the undulated topography imparts strain. Strain is known to perturb the crystal Hamiltonian through a deformation potential 19 , to a magnitude proportional to the strain, which in turn is determined by the topography of the surface. A periodic strain modulation will create a confining potential, which-if strong enough-can localize electrons and result in modulated super-lattice band whose bandwidth depends on the surface's specific geometry. Hence, patterned surfaces with specific topography can, in general, create and fine-tune flat bands in any 2D semiconductor material. To the best of our knowledge, creating either 2D or 1D flat bands in monolayer semiconductors by undulation has not been discussed before.
We illustrate this idea by theoretically investigating the electronic properties of monolayer hexagonal boron nitride (hBN) deformed by a bi-periodic sinusoidally modulated topography. Interestingly, strained hBN attains both pseudo-electric field E P (by virtue of deformation potential) and also (having a honeycomb lattice, like graphene, with two inequivalent basis atoms) pseudo magnetic field B P 20 . We find that for bi-sinusoidal deformation, E P and B P have very different spatial dependence. A combination of these leads to anisotropic confinement and creates one-dimensional (1D) flat bands, whose bandwidth can be varied by the surface topography. The 1D nature of these states can be exploited to probe the exciting physics of one-dimensional quantum systems, which have been predicted to exhibit interesting effects 21 such as Luttinger liquid behavior, charge and spin density waves, Peierls instability, and deviation from Fermi-liquid theory. The origin and nature of these 1D states are different from the 2D flat bands observed in TBG and twisted vdW heterostructures, opening an exciting realm of exploring many-body effects in 1D quantum systems in a clean and controllable manner.

Results
A 2D material conformed to a curved non-Euclidean surface undergoes a locally in-plain strain, which can be evaluated at a continuum level (for relatively smooth topography) by solving the second Föppl-von Kármán (FvK) equation 22 Here χ, Y, and f, are the Airy stress function, Young's modulus, and the surface shape function, respectively. We consider bi-periodic sinusoidally modulated off-plane shape f(x,y) = hsinαx·sinαy, with α = 2π/L, akin to egg-cart, L and h are the undulation period, and amplitude, respectively. The displacement (u x , u y ) and imparted strain (u ij ) fields were solved analytically from the FvK equation for the sinusoidal surface (for details see Methods section and Supplementary Note 1). The maximum tension or compression for this geometry depends on the aspect ratio A = h/L, as ε = |u ii,max | = π 2 A 2 (1 − ν)/2, where ν = 0.31 is the Poisson's ratio for hBN and u ii = (u xx + u yy )/2. A relaxed atomic structure of a sinusoidal hBN can be further constructed from the solved displacements u x and u y , as in Fig. 1b, for A = 0.079, where color shows u ii (for details see Supplementary Note 3). As expected, the hilltops and valleybottoms are stretched while the saddle areas are compressed, to an amplitude ε~2.12%. We note here that, generally, the material strain depends on both the surface shape f(x,y) and boundary conditions (in case of growth, also on the chemical potential μ hBN , controlled by the growth conditions). Here we allow full relaxation to minimize elastic energy while accommodating the substrate topography, that is no forces at (remote) layer's perimeter and no friction to the substrate. This corresponds to "stamping" the 2D material onto the frictionless matrix, when the layer contracts laterally, with non-negligible displacements u x , u y (see Supplementary Note 1). Analogous to graphene [23][24][25][26] , the strain in hBN generates pseudo-electric and magnetic fields 20,23 , significantly perturbing the crystal Hamiltonian. The low energy effective Hamiltonian in strained hBN in the vicinity of the K points is given by where v F = 3|t|a/2, |t| is the nearest-neighbor (NN) hopping amplitude, and a is the interatomic distance. σ (τ) =(τσ x ,σ y ,σ z ) are defined in terms of the three Pauli matrices, and τ = +1 (−1) for K (K'). k = (k x ,k y ,Δ), where ħv F Δ is the difference in sublattice potential between B and N atoms, and k x,y is the electron crystal momentum measured relative to K or K'. A p is the pseudo-vector potential caused by shear, where β 0 = (a/t)∂t/∂a = −3.3. φ p is the pseudo-electric potential (PEP) arising due to the hydrostatic component of strain, φ p = −gu ii , where g ≈ 3.66 V 27 . Accordingly, these potentials generate pseudo-electric field (PEF) E p = −∇φ p and pseudo-magnetic field and z the unit z-vector. One can already recognize that these additional pseudo fields in the Hamiltonian, arising due to strain, act as a perturbing confinement potential. The strain fields obtained for sinusoidal surfaces allow us to derive the analytical expressions for pseudo electric and magnetic fields (for details see Supplementary Note 2).
Figure 1c, d shows the PMF and PEP for A = 0.079, and L = 6.35 nm. The spatial dependence for both fields is different and surprisingly, PMF depends only on y (Fig. 1c). It is known 28 that periodic magnetic fields can lead to confinement and create localized electronic states. Similarly, we expect that for sinusoidally modulated hBN, a combination of both PEP and PMF will create flat bands. Sections of these fields along y-direction, at x = const are plotted in the right panels of Fig. 1c, d. The periodic PMF has an amplitude of B p,max~4 20T, which corresponds to confinement energy of~2μ B B p,max = 49 meV, while the periodic PEP corresponds to confinement energy of~76 meV. We will show that different spatial dependence of PMF and PEP leads to anisotropic confinement and results in the interesting electronic nature of the flat bands. We next calculate the electronic bands of our topographicallystrained hBN, using density-functional based tight-binding (DFTB) theory with a local orbital basis 29 . DFTB has been successfully applied to study various forms of hBN 15,30 , for which DFT calculations are intractable (see Supplementary Note 4 for details).
Monolayer hBN honeycomb lattice is akin to graphene, yet the different basis atoms break the sub-lattice A-B symmetry, and an energy gap opens, making hBN an insulator. The undeformed monolayer hBN shows a band gap of~3.55 eV. Figure 2a shows the band structure under bi-sinusoidal strain, with A = 0.079, ε = 2.12%, and L = 6.35 nm. The Brillouin zone is defined based on the shape function. One can see additional bands appearing in the gap, looking like defect states which might arise due to electronic confinement. We find that the bandwidth (W) of the lowest unoccupied states (shown in red) is W = 39 meV, which is very small, and it is a flat band; in comparison, the effective W of pristine BN bands corresponding to nearest neighbor hopping t~2.16 eV 20 is W~4t = 8.6 eV, which is much larger than the W of the modulated flat bands. These flat bands are well separated by >100 meV from the other states at higher energies. Interestingly, the bands are dispersive along Г-X and almost non-dispersive along R-X, which corresponds to k x and k y directions, respectively. This makes these flat bands one-dimensional and very different from those seen in TBG and other twisted vdW heterostructures. The band decomposed charge density |Ψ nk | 2 in Fig. 2b corresponds to one of the eigenstates at Г point. The electronic states are delocalized along the x-but are completely localized along the y-direction, confirming these flat bands' 1D nature (see Supplementary Fig. 6 for the charge density corresponding to the whole flat band, which is similar to Fig. 2b).
The flat bands are composed of 8 electronic states and are localized mainly on the four extremes of the sinusoidal modulation in Fig. 2b. Figure 2c shows the enlarged view of the flat bands plotted along X-Г-Y. We find that the flat bands dispersion and charge modulation can be described by a simple "coarse-grained" where t x and t y are hopping amplitudes along x-, and y-direction, respectively (for details see Supplementary Note 5). The 8 bands arise from the states localized on the 2 maxima and 2 minima, and each of them being doubly occupied. The maxima and minima act as artificial "quantum dots". Fitting this TB model to DFTB results gives |t x | = 8.9 meV (~W/4) and |t y | = 0.9 meV. The ratio of hopping along x-, and y-direction is |t y |/|t x | = 0.1, again a manifestation of one-dimensionality of the electronic states. This is quite surprising at first, because the strain pattern appears to be isotropic along the x and y-direction (Fig. 1b).
To gain microscopic insights into the reasons behind these flat bands' 1D nature, we calculated the electrostatic potential along x (Fig. 3a) or y (Fig. 3b) while averaged along the other two perpendicular directions. The sharp features are due to the approximations used to evaluate the diverging potential near atomic sites (for details see Supplementary Note 4). The potential rapid oscillations are due to periodic atomic sites, but a long-range modulation can also be seen. Gaussian averaging extracts the long-range modulation (V conf. , red lines in Fig. 3a, b), which is very different along x-and y-directions: eV conf. along y-direction has a depth of~500 meV, larger than mere~9 meV along x. This anisotropic confinement is expected because of the different spatial dependence of PEP and PMF, Fig. 1c, d. The smaller V conf. along x-direction is mainly due to contributions from PEP only, while the larger V conf. along the y-direction is contributed by both PEP and PMF. This signifies that the long-range potential modulations are due to pseudo-electric and magnetic fields, providing the anisotropic confinement needed to maintain one-dimensional flat bands. Additionally, the larger confinement energy along ydirection results in lower hopping amplitude t y in our coarsegrained model.
We find that the width of these flat bands can be tuned by changing either the aspect ratio of the topography, to alter the strain, or the period L of the undulation. Figure 4 shows the variation of W as a function of L, and ε. W goes below 10 meV for L > 8 nm. Moreover, W is found to depend exponentially on both L, and ε, a good fit W ∝ exp(−0.56*L-0.34*ε) is obtained (blue surface in Fig. 4); this is expected since W effectively corresponds to the hopping amplitude between the coarse-grained sites, defined by hybridization/overlap of wave function between them. Since the strength of hybridization decreases exponentially with distance L, W is found to show the same dependence. Additionally, the exponential dependence of W on ε can also be understood from the well-known dependence of hopping energy (t) with strain 20 , t(a) = t 0 exp(−|β 0 | (a/a 0 −1)), is the hopping amplitude at the effective bond length a, and (a/a 0 −1) is its strain. We point out that flat bands appear even at strain as low as 1.5%, and these bands are well separated by >70 meV from the other unoccupied bands above, which makes them accessible to experiments. Moreover, the quasi-1D nature of these flat bands remains robust against local strain imperfections (for details see Supplementary Note 2d) as well as small misorientation between hBN and the substrate (for details see Supplementary Note 2c).

Discussion
Interesting strongly correlated physics in one-dimension 31 is expected when the ratio of on-site Coulomb interaction U (responsible for electronic correlation) and the hopping amplitude t is large, U/ 4t = U/W > 1. Since U depends inversely on length L, U∝1/L 13 , and t∝exp(-|β 0 | L) 20 , the condition U/W > 1 should be easily achievable for reasonable L. We estimate U for our systems as U = e 2 /2πκL 13 , where L is the length of periodic modulation, and κ = 4.73 32 is the effective dielectric constant of hBN. The red surface in Fig. 4 shows values of U. U/W can be enhanced by either increasing the aspect ratio (strain) and/or increasing the periodic length (Fig. 4). The region in (L, ε) with U < W is shown by the gray shaded area in Fig. 4. To achieve strongly correlated phases, where U > W, topography with L and ε values lying outside of the gray area in Fig. 4 must be chosen. Depending on the ratio of U/W and band filling, one can expect 31 different phases such as Mott insulator (MI), Luttinger liquid, bond ordered wave (BOW), and band insulator (BI). E.g., at small band filling, gradually changing U/W from 0 to 2 can change the electronic phases in order BI → BOW → MI. Accordingly, one should expect that changing the periodic topography will provide a unique control to drive the system to different strongly correlated regimes exhibiting interesting physics. To realize the strongly correlated physics, the flat bands (Fig. 2a) must be partially filled, perhaps by electrostatic doping, as routinely done for 2D materials, including twisted bilayer TMDs 4 , and graphene 2 . Long-range ordered quantum phases in 1D tend to get destroyed due to thermal fluctuations at finite temperature 21 , hence, isolated 1D systems are not ideal to realize 1D physics. Importantly, our predicted system with parallel 1D states (resulting from the quasi  1D flat bands in a 2D material, Fig. 2b) will suppress such fluctuations, due to the finite coupling between them 21 , and will strengthen the physical effects in 1D, which will help to achieve interesting physics in 1D at finite temperatures. To realize our predictions in experiments, hBN may be overlaid or stamped (or possibly grown directly) on patterned substrate, for instance SiO 2 , having a band gap of~8.9 eV 33 , larger than hBN---so that there are no unwanted hybridization between hBNs electronic states and the substrate. Fabricating bi-periodic sinusoidal modulation on silicon is challenging but has already been attempted 34 . Delamination from the substrate, if the desired strain level is high may be a concern; by comparing the surface pressure due to substrate with the maximal adhesion forces of hBN to SiO 2 with γ ad~1 3 meV/Å 2 (see Supplementary Note 6), we estimate that a strain up to ε~2.75% at A~0.09 (see Supplementary Note 3b), should be sustainable in experiments.
In summary, we have shown that deforming a 2D semiconducting monolayer with a particular topography having nonzero Gaussian curvature, can be used as a unique and straightforward way to create flat bands and drive the system into different strongly correlated electronic regimes. Topographical modulation can be created by electron-beam lithography and does not require accurate fine-tuning of the twist angle and overcomes several challenges of twisted systems. For hBN as an example, we show that bi-periodic sinusoidal modulation generates pseudo-electric and magnetic fields, creating anisotropic electronic confinement and one-dimensional flat bands. Our proposed way to create 1D flat bands should be applicable to a variety of 2D systems like hBN with massive Dirac fermions. These flat bands are different from those observed with twisted bilayer graphene and other vdW heterostructures. These bands' one-dimensional nature will pave the route to study the exciting physics of strongly correlated 1D systems, thereby going beyond what's achievable with twisted materials. In addition to creating flat bands, substrate engineering can be used to realize intriguing electronic behaviors such as recently demonstrated electron optics and the valley Hall effect in undulated graphene 35 .

Methods
The continuum displacement and strain fields from FvK equations. For a 2D material over a gently varying topography defined as f(x,y), its structural relaxation, on a continuum level, is governed by the Föppl-von Kármán (FvK) equation 22 : where subscripts denote partial derivatives. Here χ(x, y) is the Airy stress function and Y the 2D Young's modulus. In this work we consider a sinusoidal topography f(x,y) = hsinαx.sinβy, where h defines the height, and α = 2π/L x , β = 2π/L y defines the lateral periodicity L x and L y . Plugging f(x,y) into Eq. (4) we will find And, integrating twice, we obtain the Airy function as Integration constants are set to zero to ensure the lowest elastic energy. With the Airy function, the components of the strain tensor can be obtained as u ij = (1/Y) (ε ik ε jl -νδ ik δ jl )∂ k ∂ l χ: u xx ¼ ðχ yy À νχ xx Þ=Y ¼ Àðh 2 =8Þðα 2 cos 2βy À νβ 2 cos 2αxÞ ð7Þ u yy ¼ ðχ xx À νχ yy Þ=Y ¼ Àðh 2 =8Þðβ 2 cos 2αx À να 2 cos 2βyÞ ð8Þ The continuum pseudo-electromagnetic fields from sinusoidal modulation.
Building the atomic structure. The relaxed 2D materials crystal geometry on curved surfaces was constructed using the displacement fields of the respective strain fields. This requires (1) determining the amount of material in the periodic box and (2) find the correct displacement for each atom. From the definition of the strain u ij = (1/2)(∂ i u j + ∂ j u i + ∂ i f∂ j f) we can integrate the components u xx = ∂ x u x + (f x ) 2 /2 and u yy = ∂ y u y + (f y ) 2 /2 and find the displacement fields u x , u y u x ¼ Àðh 2 α 2 =8Þx þ ðh 2 =16αÞðνβ 2 α 2 þ 2 cos 2βyÞ sin 2αx ð10aÞ u y ¼ Àðh 2 β 2 =8Þy þ ðh 2 =16βÞðνα 2 β 2 þ β 2 cos 2αxÞ sin 2βy ð10bÞ which connects the deformed coordinates (x, y) to the reference coordinates (X, Y) = (x − u x , y − u y ). The first linear term in Eq. 10 determines the overall amount of lateral contraction, and the second term is an oscillating term causing periodic stretching/compression patterns. Hence we can write (L x -L x0 )/L x = −(h 2 α 2 /8) and (L y -L y0 )/L y = −(h 2 β 2 / 8), with L x and L y the periodicity of the sinusoid and L x0 , L y0 the periodicity of the original, reference flake. Plugging in the above definitions α = 2π/L x , β = 2π/L y we have (taking i = x or y) L i 2 − L i0 L i + (h 2 π 2 /2) = 0 The atomic structure can therefore be constructed as the following: (1) Create a hBN sample with dimensions L x0 , L y0 (2) Choose desired height h for sinusoid (3) Solve for the optimal L x , L y , or equivalently α = 2π/L x , β = 2π/L y for the sinusoid (4) Displace each atom from (X, with u x , u y according to Eq. 10. DFTB calculation. The electronic structure of the pristine and sinusoidally modulated boron nitride was calculated using the density functional based tightbinding approach implemented in DFTB + using atomic orbital basis. The self consistent charge calculation was performed using a varying k-grid of 6 × 6 × 1 − 12 × 12 × 1 depending on the size of the unit cell. The maximum angular momentum chosen for B and N atoms was p (l = 1). The pairwise B-B, N-N, and B-N Slater koster files (parameterization data) were obtained from the DFTB+ 29 repository, matsci. The electrostatic potential was estimated by taking the Mulliken-point charges and superposing the corresponding 1/r potentials as implemented in the DFTB + code. In the code, the 1/r potential is modified to remove the r = 0 divergence, and instead plots 1/√r 2 + ε 2 , where ε = 10 −4 .
Coarse-grained 8 band tight-binding model. DFT calculation of adhesion energy of SiO 2 on hBN and strain energy of hBN. A 2 × 2 unit cell of hBN and 7 layers of SiO 2 along [001] was taken to create the hBN|SiO 2 slab geometry. The SiO 2 slab's surface was passivated with hydrogens. The geometry was fully relaxed using first-principles density functional theory (DFT) implemented in VASP 36 . Ion-electron interactions were represented by allelectron projector augmented wave potentials. The generalized gradient approximation (GGA) parameterized by Perdew-Burke-Ernzerhof (PBE) 37 was used to account for the electronic exchange and correlation. A plane wave basis with a kinetic energy cut-off of 500 eV was used for wave functions expansion and a Monkhorst-Pack grid of 12 × 12 × 1 k-points was used to sample the Brillouin Zone (BZ). A vacuum of 20 Å was used along the direction perpendicular to the slab to reduce the interaction between the periodic images. The structure was relaxed until the Hellmann-Feynman forces on the atoms were <0.01 eV/Å. The DFT-D2 method of Grimme was used to include van der Waals interaction.

Data availability
The authors declare that the data supporting the findings of this study are available within the paper and its Supplementary Information